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Abstract 



,_i In this paper, we study the evolution of the localized induction approximation 

Ph (LIA), also known as vortex filament equation, 

< 

-(— > 

^ for X(s, 0) a regular planar polygon. Using algebraic techniques, supported 

I *^i by full numerical simulations, we give strong evidence that X(s, t) is also a 

polygon at any rational time; moreover, it can be fully characterized, up to 
j>. a rigid movement, by a generalized quadratic Gaufi sum. 

We also study the fractal behavior of X(0,t), relating it with the so- 
called Riemann's non-differentiable function, that, as proved by S. Jaffard, 
i/~) fits with the multifractal model of U. Frisch and G. Parisi, for fully developed 

■^ turbulence. 
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'^ 1. Introduction 

Given a curve Xq : M — > M^, we consider the geometric fiow 



Xt = cb, (1) 
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where c is the curvature and b the binormal component of the Frenet-Serret 

formulae 

^T\ / c 0\ /T^ 

n = -c r ■ n I . (2) 

.b/. \0 -r Oj Vb. 

The flow can be expressed as as 

Xf = Xs A Xss, (3) 

where A is the usual cross-product, t is the time, and s is the arc-length 
parameter. It appeared for the first time in 1906 [Ij and was rederived in 
1965 by Arms and Hama [2] as an approximation of the dynamics of a vortex 
filament under the Euler equations. This model is usually known as the 
localized induction approximation (LIA); we refer the reader to [3] and jl] 
for an analysis and discussion of its limitations. 

([I| and ([3]) are also known as the binormal equation and the vortex fil- 
ament equation, respectively. Some of their explicit solutions are the line, 
the circle, and the helix. Since the tangent vector T = X^ remains with 
constant length, we can assume that T e S^, Vt. Differentiating ([s]), we get 
the so-called Schrodinger map equation on the sphere: 

Ti = T A T,„ (4) 

which is a particular case of the Landau-Lifshitz equation for ferromagnetism 
[5]. Q can be rewritten in a more geometric way as 

T, = JD,T„ (5) 

where D is the covariant derivative, and J is the complex structure of the 
sphere. Written in this way, Q can be generalized to more general definition 
domains and images, as the hyperbolic plane H^. 

([3]) is time reversible, i.e., if X(s, t) is a solution, so is X(— s, —t). Bearing 
in mind this fact, an important property of ^ (and, hence, of Q) is that 
it has a one-parameter family of regular self-similar solutions that develop 
a corner-shaped singularity at finite time. This was proved in [6j for the 
Euclidean case, and in [7], for the hyperbolic case. Furthermore, starting 
from T. F. Buttke's work |8], there have been a couple of papers devoted to 
studying numerically the self-similar solutions of (IS]) , as in [H] , where a careful 
numerical study of those solutions for both the Euclidean and the hyperbolic 



cases was done. In [H], on the one hand, the authors reproduced numerically 
the formation of the corner-shaped singularity, and, on the other hand, they 
started with a corner-shaped initial datum, recovering numerically the self- 
similar solutions; in all cases, the correct choice of boundary conditions was 
shown to be vital. Furthermore, they gave numerical evidence that fractal- 
ity phenomena appeared in (|4]), if, for instance, fixed boundary conditions 
were imposed on T. Remark that the relationship between fractals and the 
Schrodinger map is not new; indeed, in [10], an aortic valve model was pro- 
posed, in order to study the apparent fractal character of the valve's fiber 
architecture. Instead of X(s,t), the authors wrote X(-u,f), where v, which 
corresponds to our time, is such that the curves v = constant are the fibers. 
After imposing fixed boundary conditions at u = ±uo, they showed that the 
curves u = constant have a fractal character. Besides, the fractal dimension 
of those curves was calculated numerically in |llj . 

Even if the solutions of (tsl) for an initial datum with a corner are well 
understood, nothing had been done for more general initial data with several 
corners, in particular, polygons. Nevertheless, in a recent submitted paper 
|12j by R. Jerrard and D. Smets, they propose a global existence theorem 
that allows to consider such types of initial data. Moreover, they simulate 
numerically the evolution of the unit square at different times, suggesting 
that the solution could become again polygonal; indeed, at a certain time, 
the square seems to reappear, but rotated 7r/4 with respect to the initial 
one. In the following pages, we will show not only that these observations 
are correct, but that, given a regular planar polygon as initial datum, we 
have a completely describable polygon at any rational time. 

The structure of this paper is as follows. In Section [2| we apply the 
Hasimoto transformation to Q, but considering a variation of the Frenet- 
Serret formulae (|2]). This relates (tsl) and ^ with the nonlinear Schrodinger 
(NLS) equation: 

ij, = ii:^^ + i(^{\ij\^ + Ait))\i:, A{t)eR, (6) 

where A{t) G M can be written in terms of c(0, t), r(0, t) and their derivatives. 
The fact that ([3])-(|4]) and the NLS equation are related has far-reaching 
consequences. 

In Section [3| we gather all the theoretical arguments that support our 
numerical experiments. We start recalling some elemental geometric sym- 
metries that regular polygons have and that are preserved by LIA. On the 



one hand, these symmetries greatly simphfy the numerical implementation 
and, on the other hand, they play a fundamental role in computing some 
important quantities that govern the dynamics of regular polygons. Then, 
we pay special attention to a group of symmetries that, as far as we know, 
can only be visualized through the use of the Hasimoto transformation and 
the NLS equation (|6]). These are the so-called Galilean transformations (33), 
that leave invariant the set of solutions of LIA. 

We will proceed as in [13] , where the solution of NLS that has the delta 
function as an initial condition is characterized by the fact that it leaves 
invariant all the Galilean transformations. This property determines the 
solution, except for one function of time that is easily computed by solving a 
simple ODE and choosing A{t) in (p| appropriately. This argument heavily 
relies on assuming that uniqueness holds for such a singular initial datum. 
This is in indeed a quite delicate question that in the case of just one delta 
(i.e. a curve with just one corner) has needed some work to be eventually 
answered in the positive [HI [13 [IS] • However, in this paper, we are interested 
in just proving the existence of a solution starting with a polygon. Therefore, 
we can assume uniqueness to conclude that the expected solution has as many 
symmetries as possible and simplify in this way the analysis. 

In our case, we take a planar regular polygon of M sides as initial datum. 
In terms of NLS, this amounts to considering periodic boundary conditions 
and an initial datum which is given by an infinite sum of delta functions with 
appropriate weights. Therefore, just a subgroup of the group of the Galilean 
transformations has to be considered. Then, as in |jl3j, the solution of (|6| is 
determined except, again, for a function that will depend on time. But, in 
this case, it is very delicate how to determine this function and does not seem 
to be easy. Moreover, a very simple analysis makes clear that it cannot be 
expected to make any sense to the NLS equation ([6]). In fact, the arguments 
in [16j prove that, in the case of an initial condition given by a curve with 
just one corner, ^ does not make sense in an interval of time that contains 
t = 0. Our approach in this paper is to integrate the Frenet-Serret system 



(18) with ip{s, t) = a + i/3 the distribution given by 

oo 

^{s,t) = V'(0,t) Y, e-^(*^'=)'*+'*^^^ (7) 



We do it just for times that are rational multiples of 27r/M^, that is to say, of 
the form tpq = {2n / M"^) {p / q) , because in that case the corresponding curve 



is a skew polygon with Mq sides, for q odd, and with Mq/2 sides, for q even, 
determined by a generahzed quadratic Gauss sum: 



c-l 

;2 



G(a,fe,c) = ^e2'^'('^''+'')/^ a,5eZ,cGZ-{0}. (8) 

z=o 

The function of time ip{0,t) is then obtained by imposing the condition that 
the polygon has to be closed. As a consequence, the tangent vector is fixed, 
except for a rotation that can depend on time. This rotation is determined 
by using the extra symmetries that are available for regular polygons and 
that we have already mentioned. The next step is to integrate T, in order 
to obtain X at rational times. Using again the symmetries, we are able to 
calculate all the necessary quantities, except for some possible translation in 
the z direction that can be easily computed numerically. Nevertheless, we 
do not see a simple way of obtaining this last value using just theoretical 
arguments. 

At this point it is important to recall the invariances of the so-called 
Jacobi theta function, which is closely related to rt7|: 

oo 

e{s,t) = Y, e-^'^^'^^^''"' . (9) 

fc=— oo 

This function is precisely the solution of the linear Schrodinger equation 

i't = ^i^ss, (10) 

with a 1-periodic delta as the initial condition. Therefore the Galilean trans- 
formations still hold in this case; moreover, they also imply the extra sym- 
metry 

e{s,t) = e^-/'^-^e('-,=p), t>o, (11) 



Vi \t' t 
which, together with the periodicity property, 

e{s,t + 2) = e{s,t), (12) 

generates the so-called unimodular subgroup. In fact, in [T71 [TS], all these 
symmetries are used to prove the existence of a continuum range of exponents 



(usually known as the spectrum of singularities) for the Holder regularity of 

'^W = E|^' ^^[0'2]. (13) 

fe=i 

whose real part is precisely Riemann's non-differentiable function: 

k=l 

More concretely, in [THj it is proved that Riemann's non-differentiable func- 
tion is a multifractal. In other words, it is proved that the set of times t 
that have the same Holder exponent is a fractal with a dimension depend- 
ing on the Holder exponent, in such a way that the conjecture stated by U. 
Frisch and G. Parisi in [12] is fulfilled (see also [2D], for more details at this 
respect and the connection of this question with fully developed turbulence 
and intermittency) . 

From all this, it is very natural to ask if something similar happens with 
the evolution of a regular polygon by LIA. There are different ways of under- 
standing the question. One is whether, for a fixed time, the corresponding 
curve is a fractal or even a multifractal. From what we have said above, this 
cannot happen for rational times, because, in this case, it is a skew polygon. 
But precisely these examples do suggest that quite likely the situation for 
irrational times is completely different; and in fact it is, as it is explained 



in Section 5.2, At this respect, we have to call the attention of [21], where 
upper bounds for the regularity of (|9]) measured in Besov spaces of negative 
indices are obtained. 

On the other hand, we could also ask if something similar to the properties 
of Riemann's non-differentiable function also holds in our case. It is very easy 



to find the analogous to (13) in out setting. In fact, for a given M, if we 
write X(0,t) = (Xi(0, t),X2(0, t),X3(0, t)), bearing in mind the symmetries 
of the problem, we conclude that X(0, t) is a planar curve, so we identify the 
plane where it lives with C and define 

z{t) = -||(Xi(0,t),X2(0,t))|| +zX3(0,t). (15) 

This curve, or, more precisely, z{t) — CMt, for a certain cm > that depends 



on M, can be seen as nonlinear versions of (13). Is z{t) a multifractal? If the 



answer is positive, what is its spectrum of singularities? We consider that 



these two questions are rather challenging from the numerical and analytical 
point of view. The results of Section |4] of this paper strongly suggest that 
the answer to the first question is positive. In fact, we compute z{t) — 
CMt for different M and measure the error of the difference between it and 
an appropriately scaled and rotated version of 0(t) in the L°° norm. The 
convergence rate is rather strong. 

We have also made some numerical experiments that prove that the 
Holder exponent of z{t) is 1/2 for rational times tpq = (27r/M^)(p/g). Nev- 
ertheless, how the constants depend of the denominator q, something which 
is a fundamental ingredient in the arguments in [T7] and in [TS], is unclear. 
This question deserves a much more detailed analysis that we plan to make 
in a forthcoming paper. 

In Section 111 we propose a numerical method for ^ and (El). We explain 
carefully how to take advantage of the symmetries of X and T for the types 
of solutions considered, which leads to a dramatically improvement in the 
computational cost of the algorithm. 

In Section [5} we perform numerical experiments, simulating the evolution 
of ([3]) and ^ for different numbers of sides M in the initial polygon. We 
compare the numerical results with the algebraic values from Section [3} ob- 
taining a complete agreement between both (totally different) approaches. 
We also analyze carefully the evolution of X, for s = 0, giving numerical ev- 



idence that and z{t) in (13) and (15) are intimately connected and, hence. 



that z{t) is also a multifractal. Finally, we make some comments on the 
structure of T(s, t) at a fixed time tpq with g :» 1. 

In Section [6} we offer the main conclusions, as well as some open questions 
that we postpone for the future. 



In order to conclude this paper, we offer in Appendix A a detailed study 
of those aspects of the generalized quadratic GauB sums that are of relevance 
for our work. 

2. The Hasimoto transformation 

A central point of this paper is the natural connection between (pi)- (111) and 
the nonlinear Schrodinger (NLS) equation. Indeed, applying the Hasimoto 
transformation l22l: 



ip{s,t) = c(s,t)exp ( i / T{s',t)ds' ) , (16) 



ip satisfies the equation 

ijjt = iipss + i 



+ Ait)) up 



(17) 



where A{t) is a certain real constant that depends on time. Nevertheless, 
for our purposes it is not convenient to work with the torsion. Instead, we 
consider another version of the Frenet-Serret trihedron; it is easy to check 
that all its possible generalizations have the form 




a I3\ T' 
-a 7 ■ I ei 
-/3 -7 0/ Ve2, 



(18) 



for some vectors ei and 62 that form an orthonormal base with T. Moreover, 
we can make, without loss of generality, one of the coefficients a, /3, or 7 equal 
zero. If we make /3 = 0, denoting a = c and 7 = r, we recover the habitual 
trihedron. On the other hand, in this paper, we choose 7 = 0, in order to 
avoid working with the torsion. In that case, the Hasimoto transformation 
takes the form 

ij = a + i(3. (19) 

It is quite straightforward to check that this new definition of ip also satisfies 



(17). Indeed, defining (see [22] 



N = ei +«e2, 
the generalized trihedron becomes 



^(^N + ^N) 



Bearing in mind that T A N = — zN, Q gets transformed into 

"1 



(20) 



(21) 



Tt = TA 
1 



-(^N + ^N) 



T A (^,N + ^N, + ^,N + i)'N, 
(^,N-^,N). 



(22) 
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We want to express N^ as N^ = aT + 6N + cN: 

(23) 



N • T = ^ Nt ■ T = -N ■ Tt = - 


-itjj 


N-N = O^NfN = 




N-N = 2^NfN + N-Ni = 0. 





From the last equation, it follows that the real part of Nj ■ N is equal to zero, 
so Nj ■ N = iR, for some R ^M.. Putting all together, 

Ni = -i^,T + iRN. (24) 

Differentiating this last expression with respect to s, 

Nt, = (-2^,T + tRN), 

= -i^pssT - iipsTs + iRsN + iRNs 

= -(#,, + iRij)T - '-ijiJsN + zRsN - ^V^V^.N. (25) 

On the other hand 

N,i = -(V'T)i 

= -^tT - ^Tt 

= -^tT + ^V'V'sN - '-^^s^. (26) 



Equating the coefficients of T and N in (25) and (26): 



1 - 1 - 1 (27) 

Rs = -^^s + ^^^s ^R= -m^ + A{t)), A{t) G M. 



This concludes the proof that ■?/' = a + i/3 satisfies (17). 



3. A solution of Xt = Xg A X^s for a regular polygon 

The aim of this paper is trying to understand the evolution of (tsl) for 
polygonal initial data. In this paper, we limit ourselves to studying the 
simplest case, i.e., that of a regular planar polygon of M sides. Since ^ and 
^ are invariant by rotations, we can assume without loss of generality that 
X(s, 0), and hence T(s, 0), live in the plane OXY, i.e., their third component 



is zero. Identifying the plane OXY with C, and assuming without loss of 
generality a total length of 27r, X(s, 0) is the polygon parameterized by arc- 
length whose M vertices sq, . . . , sm_i, located at Sk = 2nk/M, are 

^'--■°'- Msin(./M) - P«) 

and X(s,0), for Sk < s < Sk+i, is in the segment that joints X(sfc,0) and 
X(sfc+i,0). Observe that it is straightforward to extend X(s,0) periodi- 
cally to the whole real with period 27r. Then, since X(sfc+i,0) — X(sfc,0) = 
(27r/M)e^'^*'^/*^, the corresponding tangent vector is the periodic function 
with period 27r such that 

T{s, 0) = e^"^'/^, for s^ < s < Sk+i. (29) 

Remark that X(s, 0) is continuous, while T(s, 0) is only piecewise continuous. 
In fact, X(s, 0) can be regarded as a curve whose curvature, periodic with 
period 2'k/M, is zero everywhere, except at the vertices, where it equals 
infinity, i.e., it can be expressed as a sum of Dirac deltas: 



oo 



<-') = M E^(^-¥)' (30) 

A;=— oo 

where the constant 2tt /M has be chosen in order that the integral of the 
curvature over an interval of length 2tt is equal to 27r: 



27r- 271- l"i-^' 

n{s)ds = — - 



'M-l 

27rfc ^ 
M > 
k=0 



E^(»- 



ds = 27r. (31) 



At this point, it is vital to understand the role played by symmetries in ([s]) 
and ^ for this initial datum. Both equations are invariant by rotations, 
i.e, given a rotation matrix R, if X = (Xi,X2,X3)^ and T = (Ti,T2,T3)^ 
are, respectively, solutions of them, so are R ■ X and R ■ T. Therefore, 
if R ■ X(s,0) = X(s,0) and R ■ T(s,0) = T(s,0), then, if the solution is 
unique, R ■ X(s,t) = X(s,t) and R • T(s,t) = T(s,t), Vt. In particular, 



since X(s, 0) and T(s,0) defined in (28) and (29) are invariant to rotations 



of angle 2Trk/M, around the 2;-axis, VA; G Z, we conclude that also X(s, t) 



10 



and T(s, t) are invariant to that kind of rotations, Vt. More precisely, writing 
X = (Xi,X2,X3), T = (Ti,T2,T3), then, VA; G N, 

Xi(s + ^,t) +2X2(5 + ^,t) = e2-'=/^^(Xi(s,t) +2X2(5, t)), 

X3(5 + ^,t)=X3(5,t), 

Ti(5 + ^,t) + zT2(5 + ^,t) = e'-^'/^\T^{s,t)+tX,{T,t)), 

consequently, for a given t, X(s + 27ik/M,t) lay in the same orthogonal 
plane to the 2;-axis, VA; G N. This fact has very important implications for 
implementing efficient numerical schemes, as we will see in Section |4| 

Besides invariance by rotations, (|3| is also mirror invariant. Indeed, it 
is straightforward to check that, if, for instance, X is a solution of (pi), so 
is X(s,t) = (-Xi(-s,t),X2(-s,t),X3(-s,t))^, etc. Therefore, if X(s, 0) = 
X(s,0), then that symmetry will be again preserved during the evolution, 
i.e, X(s,t) = X(s,t), Vt, or, in other words, X(s,t) and X(— s,t) will be 
symmetric with respect to the plain containing the z-axis and the y-axis, 
Vt. In our case, we can see immediately that a regular polygon of M sides 
like (28) has 2M of those mirror symmetries. For instance, in Figure [11 we 



have plotted a pentagon (blue), together with the projections over z = of 
its ten symmetry planes (red). Those symmetry planes are such that they 
contain the z axis, and one of the points X(7r/c/5, 0), A; = 0, . . . , 9, which are 
precisely the vertices of the pentagon and the middle points of the sides of the 
pentagon. Therefore, it follows that X(7rfc/5, t) will remain in the plane that 
contains X(7r/c/5, 0) and the 2;-axis, Vt. An extremely important corollary 
that will be used later is that X(27r/M, t) — X(0,t) is a positive multiple of 
the vector (1, 0, 0)^. 



Another central point of this paper is the fact that the NLS equation (17) 



is invariant by the so-called Galilean transformations, i.e., if t/^ is a solution 



of (17), so is 

4(s, t) = e'^'-'^^'%l){s - 2kt, t), Vfc, t G R. (33) 

Therefore, if we chose an initial datum such that ipk{s, 0) = ip{s, 0), VA; G M, 
i.e., such that ip{s,0) = e^'^^ip{s,0), VA; G M, then, if the solution is unique, 
7/>(s, t) = e'^'-'^''^^{s - 2kt, t), for VA;, t G M. 

In the case of a planar initial datum X(s, 0) of ([3]), the torsion is always 
equal to zero and, V'(s,0) is by definition the curvature of X(s,0). Hence, 

11 
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Figure 1: The ten mirror symmetries for a pentagonal initial datum 



the '?/'(s,0) corresponding to a regular polygon of A^ sides is given by (30), 
i.e., 



V'(.,0) = 9E*( 



s — 



2TTk \ 

M >■ 



(34) 



fc=— oo 



Like (30), ^'(■^^O) is periodic of period 2tx/M and, since (17) is invariant by 



space translations, ip{s,t) is also periodic of that period, Vt. Moreover, it 
satisfies '?/'(s,0) = e**^^**'?/'(s, 0), V/c G Z. Therefore, the Galilean transforma- 
tion hold, i.e., V^(s,t) = e^*-^'^"-^(*^'^)'*^(s-2MH, t), \/k e Z, Vt G M. Bearing 
in mind the previous arguments, let us calculate the Fourier coefficients of 

^(s,t): 



4^{j,t) 



-iMjt 



-iMjs 



2^ 
2^ 

2^ 



ip{s,t)ds 



.JMks-i(Mk)h^^^_2Mkt,t) 

e-iM(j-k)s^^^_2Mkt,t)ds 



ds 



-iM(j-k)(s+2Mkt) 



2vr JO 

j}^^-i{Mkft-iM{j-k){2Mkt) r-2-K/M 

"^ Jo 

^_ii^Mkft-iM{j-k){2Mkt)^(^- _ ^^^y 



4'{s,t)ds 



-'m^-f'>^P(s,t)ds 



(35) 
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This identity holds for all j and k. In particular, evaluating both sides at 
j = k: 

^(A;,t) = e-*(*''=)'*V'(0,t), (36) 



so i/j can be expressed as 



iP{s,t) = 4j{0,t) Y, e^ 



i{Mkft+iMks 



(37) 



k=—oo 



where ip{0,t) is a constant that depends on time and that has to be chosen 
in such a way that the corresponding X and T are periodic of period 27r for 
a fixed t; for instance, when t = 0, we have trivially ip{0, 0) = 1. Moreover, 



combining (34) and (37), we get the following well-known identity: 



y^ gi(Mfc)s ^ 



27r 
M 



5:*(' 



27^fc ^ 
M > 



(38) 



Observe that, if we solve the linear Schrodinger equation 



Itj 'iUxx 



for the initial data (34), we get 



(39) 



Hs,t)= J2 



-i{Mkft+i{Mk)s 



M M2- 



Ql M M^\ 



(40) 



where 9{s,t) is the well-known Jacobi theta function (|9|. (40) is the same 
formula as in (37), but with '(/'(O, t) = 1, Vt G M. Indeed, obtaining a uni- 



versal expression for ^{0,t) in the nonlinear case is not an easy problem. 
Nevertheless, after observing that e"**^*^^-* *+*(^^'=)* is periodic in time with 
period 27r/M^, we have been able to calculate ilj{0,t) by algebraic means 
for rational multiples of 2-7? /M^, as is shown in (87). Moreover, that value 



appears to be unique. Therefore, we give evidence that also V'(s,t) is peri- 
odic in time with period 27r/M^, strongly suggesting that ^ is also periodic 
in time with the same period. These considerations are fully supported by 
numerical simulation. 
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3.1. ip{s,t) for rational multiples of t = 27r/M^ 

Let us evaluate (37) at t = tpg = (27r/M^)(p/g), where p eIj, q E'H, and 



we can suppose without loss of generality that gcd(p, g) = 1: 

oo 
^(s,tp,)=V>(0,tp,) Y, ^-^iMkf2.p/i^P,)+^Mks 

k=—oo 

oo 

fe=— oo 
q—1 oo 

1=0 fc=— oo 
9-1 

= 7/'(0,tpg)^e-2-'(^/'')''+**^'^ ^ e'*'^'=^ (41) 



fe=— oo 

oo 

^ — ^HLyy/^jh - 
fe = — OO 

q—1 oo 

'27rJ(p/(5r)(qfc+i)2+iM((jrfe+0'5 

/=0 fc=— oo 

9—1 oo 



/=0 fe=— oo 



Using the identity (38), 

g-i oo 

^ Z=0 fc=-oo 

^ 1=0 k=-oa 

r, oo 9-1 

^ fc=-oo 1=0 

oo (j-1 g-1 

y fe=-oom=0 1=0 

27r "^ ''"^ 

mu Yl E^(-P'^'^)^(^-^-w)' (42) 



Mg 

fe=— oo m=0 



where 



c-l 

;2 



G(a,6,c) = J]e2'^'('^''+''')/^ a,6 6 Z,c e Z - {0}, (43) 



=0 



denotes a generalized quadratic GauB sum. In Appendix A\ we give all the 



information about these sums relevant for this paper; in particular, we show 
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that, if gcd(p, q) = 1, then 

{y/q, if q is odd, 
A/2g, if q is even and q/2 = m mod 2, (44) 
0, if q is even and q/2 ^ m mod 2. 

Therefore, ip{s, t) has evolved from M Dirac dehas in [0, 2n) at t = 0, to Mq 
deltas in [0, 27r) at tpq, for q odd; and Mq/2 deltas in [0, 27r) at tpg, for q even. 



Furthermore, from (44) 




if q is odd, 
G{—p, m,q) = ^ i/Sge*^™, if q is even and q/2 = m mod 2, (45) 

if g is even and q/2 ^ m mod 2, 

for a certain angle 9m that depends on m (and, of course, on p and g, too). 



Introducing (45) into (42), and restricting ourselves to fc = 0, i.e., [0, jj), we 



conclude that 

' 7^,m U Erio e''-5is - ^), if g odd. 



i>{s,tpq) 



Ma 



V^(0,tp,)Ero e^'^'""^'^(^-^^^), ifg/2odd, ^4g) 



7771^(0, t„) E'rii'o' e^'^-^i^ - ^), if g/2 even. 

- V 2 

The coefficients multiplying the Dirac deltas are in general not real, except 
for t = and ti^2 = tt/M'^. Therefore, ipi^^tpq) does not correspond to a 
planar polygon, but to a skew polygon with Mq sides, for q odd; and to a 
skew polygon with Mq/2 sides, for q even. Moreover, ^^(O,^^^) has to be 
determined in such a way that the polygon is closed; this is done in the next 
section, where we show that its rather involved choice is unique, concluding 
that ip{s,t) is periodic in time, with period 27r/M^. 

Since the Dirac deltas are equally spaced at a time t = tpq, the length of 
the sides is the same. Observe that, if g/2 is odd, i.e., if g = 2 mod 4, then 
there is no vertex at s = 0. For instance, if p = 1, g = 2, i.e., t = tt/M^, then 

i'(».M = |f:*(»-2^^i7^). (47) 

fc=— oo 

Furthermore, the coefficients multiplying the Dirac deltas are now real, so 
"^("5,^1,2) is simply the curvature, and V'(0,ti,2) = 1 has been chosen in order 
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that the integral of ip{s,ti^2) over [0, 2tt) is 27r. Bearing in mind the symme- 
tries of the problem, we conclude that the corresponding polygon is a regular 
planar polygon of M sides, which lives at a certain plane z = constant, and 
which has been rotated n/M degrees around the z-axis with respect to the 
initial data, as R. Jerrard and D. Smets predicted. 

3.2. Recovering X and T from ip at t = tpq 

Given ip{s,t) = a{s,t) + i/3{s,t), recovering T, ei and 62 from ip implies 
integrating 

'T\ / a /3^ 

eA = l-a I • I ei I . (48) 

As we have seen in the previous section, at a time tpg = (27r/M^)(p/g), 
ip^s, tpq) is a sum of Mq (if q odd) or Mq/2 (if q even) equally spaced Dirac 
deltas, that corresponds to a skew polygon X.{s,tpq) of Mq or Mq/2 sides. 






ad h6\ 


/T 


-a6 





ei 


-b5 


0/ 


\e2 



To integrate (48), we have to understand the transition from one side of the 
polygon to the next one. In order to do that, we reduce ourselves, without 
loss of generality, to a certain if) formed by a single Dirac delta located at 
s = 0, i.e, ip{s) = (a + ib)6{s), so we have to integrate 



(49) 



This system corresponds to some T(s), ei(s) and 62(5) constant, except 
at s = 0, where there is a singularity. Let us suppose that, for s < 0, 
T(s) = T(O^), ei(s) = ei(0"), 62(5) = 62(0^), so we want to calculate 
T(0+), ei(0+), 62(0+), such that, for s > 0, T(s) = T(0+), ei(s) = ei(0+), 
62(5) = 62(0"*"). In this case, the problem is equivalent to that one of a 
plane curve with curvature k the delta function at one point, k,{s) = aS{s); 
therefore, it is equivalent to solving 




T\ _ / a6\ /T 
n ~ [-a6 J ■ I n 



(50) 



The correct way to understand (50) is to write it as 



^fo-i^')ds'^'^s)y = (e^'^/o ''(-')'^^'/(s))' = 0, (51) 
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with z{s) = x{s) + iy{s) the parametrization of the curve with arc- length 
parameter s. 



In order to integrate (49), it is important to notice that, ahhough T, 



ei and 62 are vectors, it can be uncoupled into three systems of ODE's 
for (Ti, 61,1,62,1), (T2, 61,2, 62,2) and (T3, 61,3, 62,3), respectively. Hence, let us 
consider 







ad 


h6\ 


(ui 


-a6 








■ U2 


M 





0/ 


\U3 



5(8) 



a b\ 


/ui 


-a 


U2 


-b 0/ 


\uz 



(52) 



Then, motivated by (51), we proceed as in the case of a plane curve: 




exp 



exp 




We can compute explicitly the matrix exponential: 

/ 

exp (A) = 



a sin(\/a?+6^) 



cos(v'a^ + V^) 
— asin(\/a^+b^) a? cos(\/a?+W)+b'^ abcos{\/a?+E^)~ab 






Vaa+foZ 



a2+62 



a2+62 



V 



-bam(^/a?+ifi) abcos{\/a^+b'^)—ab b'^ cos{\/a^T^)+a'^ 
Va2+b2 ^2q:^2 ^2^ 



(53) 



(54) 



Furthermore, this matrix can be rewritten in a much more elegant way. In- 
deed, expressing a + ib in polar form, i.e., a + ib = pe^^ , then 



exp(A) 



cos(p) sin(p) cos(S) sin(p) sin{6) 

- sin(p) cos{9) cos(p) cos2(6') + sm'^{6) [cos(p) — 1] cos(0) sin(S) 

- sin(p) sm(9) [cos(p) — 1] cos{9) sm{9) cos(p) sui'^{9) + cos'^{9) 



(55) 



I + sin(p) 
+ (1 



008(6*) sin(6l) 

^ cos(6») 

^sin(6») 

-1 





cos(p)) ( -cos2(e) 

- cos(6i) sin(6i) 



- cos(6l) sin(6») 
-sin2(0) 



(56) 



Hence, exp(A) is a rotation matrix [23], because exp(A)-^-exp(A) = I. More 
precisely, it is the rotation matrix corresponding to a rotation about an axis 
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(0, sin(6'), — cos(^)) by an angle p. exp(A) also relates {T(0 ),ei(0 ),e2(0 )} 
and{T(0+),ei(0+),e2(0+)}: 



exp(A) • e^iO-f . (57) 





Therefore, if {T(0 ),ei(0 ),e2(0 )} forms an orthonormal basis of M.^, so 
does {T(0+),ei(0+), 62(0+)}. 



Coming back to the general form of ^p, we have to integrate (52) Mq or 



Mq/2 times to obtain a closed skew, i.e., non-planar polygon with Mq or 



Mq/2 sides. But, according to (46), in [0 



27^^ 



'' Af /' 



Tm+2 
Mq 



i^{s,t„) = I Et=o\^2n.+i + tP2:m+i)S{s - ^^^), if q/2 odd, (58) 



Xti~o\(^2m + ^P2m)S{s - ^), if q/2 even, 



where 



2-ir 

2tt 



ip{0,tpq), if q is odd. 



\am+il3m\ = P= { 777? ^(0' ^p?)' if ^ is even and q/2 = m mod 2, (59) 

V 2 
0, if q is even and g/2 ^ r/i mod 2, 

so we conclude that, at any time tpg, the angle p between two adjacent 
sides is constant. Furthermore, the structure of the polygon is completely 
determined by the angles 9m, appearing in the generalized quadratic Gaussian 
sum, where am + i(3m = pe*^™ . 

Let be M^ the rotation matrix corresponding to (a^ + iPm)^- If «m + 
iPm = 0, Mm is simply an identity matrix and can be ignored. Otherwise, 



from (55) 



cos(p) sin(p)cos(6'm) sin(p) sin(e,„) 

Mm = \ - sm(p) cos{em) cos{p) cos^ (6*™) + sin^ (6*™) [cos{p) - 1] cos(em) sin(em) | . (60) 

- sin{p) sin{em) [cos(p) - 1] cos(6'm) sin(em) cos{p) siii^ (9^,) + cos2(9„)^ 

If, for instance, q is odd, bearing in mind that T, ei and 62 are piecewise 



constant, we have 






■ 27r 
•Mq 



27r - 
Mq 



\ e2( 







Mn 




/ T(|^ + f \ 



ei 



Ait 
M2_ 



\ 62 



' Ait 
^Mq 






eii 



■Mq 
' 27r - 



ei 






Ve2(CF" 






\ 62 



Mg 



61' 



A-K +\r 



Ml 



Mo 



61' 



27r - 
Mg 



\T 



V e2(|^-f / 



V e2(S-)^ / 



6l' 



47r ' 
Mq 



\T 



\ ^^iw-q~y I 



and so forth, i.e., there is a jump at s 



ll^. This is equivalent to writing 



rp/ 27rfc 



\ 



6l 



'27rfc + \T 
^Mg / 
2iTk + \T 



T(0- 



\T 



Mfc-Mfc_i-...Mi-Mo 



6i(0- 



(61) 



(2TTk^\T / 



62(0 



-\T 



where k G Z+, and M^ is periodic modulo q, i.e., M^+g = M^. This formula 
is also valid for q even, provided that we bear in mind that half of the M^ 
matrices are simply identity matrices. More precisely, for g = mod 4, then 



/ rp/Ank + \T 
I ^^ Mq I 

A-Kk + \T 



\ 



6l( 



Mq 



T(0- 



M2fc-M2fc-2-...M2-Mo- 61(0- 



(62) 



V ^r^T J 



62(0' 



V/c G Z+. In this case, the matrices with odd subscript are identity matrices, 
and 



■ Mq 
' A-Kk+2iT ' 
^ Mq 



V62 



■A-nk+2iT + \T 



Y J 



Mq 

Equivalently, if g = 2 mod 4, then 

' A-Kk-2TT + \T 



1 


rp ( ATTk+2lT ~ 

^y Mq 


T 


\ 


(A-Kk+2-K 
*^U Mq 


T 


\ 


(A-nk+2TT 
^^21 Mq 


T 


J 



( T(g-^)^ \ 



61 



^47rfc + ^T" 
-Mq 



(63) 



V e2(lfF / 



61I 



\ 62 



- Mq I 
' A-Kk-2-K + \T 
~ Mq I 



-*- \Mq 



M 



2fc+l 



M 



2k~l 



M3 ■ Ml 



62 1 



27r 

^MjL. 
' 2-K - 
-Mq 



(64) 
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V/c G Z+. In this case, the matrices with even subscript are identity matrices, 
and 

'4TTk-\T \ / rp/ 47rfc-27r + '\T \ 



/ rp/ 47rfc " 



\ 






T( 



Afq 



'47rfe" 



e2 



•47rfc" 



r J 



Mq 



'4T7k-2TT+\T 



ei 



Mg 



\e2{ 



4iTk-2 n + \T 
Mq 



(65) 



)^ / 



When q is even, (63) and (65) explain why the polygon has g/2 sides. Equiv- 



alently, it could be regarded as having q sides, but half of them being in- 
distinguishable, because when the angle between two adjacent sides is zero, 
they merge into a single side. 



In order that the polygon is closed, we have to choose p in (59) in such a 
way that T, ei and e2 are periodic, i.e.. 



T(27r 



-\T 



T(0 



-\T 



ei(27r 



-\T 



ei(0 



-\T 



-\T 



e2(27r 
which is equivalent to imposing that 



62(0 



-\r 



(66) 



M 



Mq-l 



■M 



Mq-2 



. . . • Ml ■ Mo = I. 



Let us define 



M = M 



g-l 



M 



q-2 



■Ml -Mq. 



(67) 



(68) 



From (67), M is an M-th root of the identity matrix. Moreover, it is also a 



rotation matrix that induces a rotation of 2n/M degrees around a certain ro- 
tation axis. Therefore, we have to choose p in order that any of the following 
properties is satisfied: 



Tr(M) = l + 2cos(|f; 



A(M) = {l,e 



2TTi/M -2■Ki/M^ 

1^ J 5 



% 



(69) 
(70) 



were Tr(M) and A(M) denote the trace and the spectrum of M, respectively. 
We work with the trace, because it is algebraically easier. In order to un- 
derstand its structure, we analyze a few cases. Let us start with q odd. The 
angles corresponding to the generalized quadratic GauB sums for the simplest 
case, g = 3, are 



P 
P 



1 : 

2 : 



^0 



2' 



01 



6' 



Oi 



02 
O2 



(71) 
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Bearing in mind that M = M2 ■ Mi ■ Mo, a symbolic manipulator yields for 
both p = 1 and p = 2 the following expression: 



Tr(M) 



;i + cos(p))= 



1 + 2 cos 






SO the only possible real value for cos(p) is 

cos(p) = (4 + 4cos(|f))i/2-l = 2cos=^/^( 



M' 



(72) 



(73) 



Notice that we arrive at the same conclusion if we impose that C = ^ 
an eigenvalue of M, i.e., if det(M — ^I) = 0: 

= det(M - CI) 



2Tri/M ■ 



IS 



c^-c 



(cos^(p) + 3cos^(p) + 3cos(p) - 2C - 3 - 2C~^) 
((l + cos(p))3-2C-4-2r^). 



Hence, 



cos(p) = (2C + 2C-^+ 4)1/3 _i 

= (2e2-*/^ + 26"=^"*/^ + 4)1/=^ - 1 
= (4 + 4cos(||))i/3-l 



2cos^/'(§; 



1. 



Let us consider now the case q = 5. Its corresponding angles are: 



= 1 


Oo = 0, 01 = 


= 2 


00 = -vr, 01 


= 3 


00 = -vr, 01 


= 4 


00 = 0, 01 = 



27r 
' 5 ' 



02 = 
5, 02 
, 02 = 



2n 
5 







5' "^ 5 

2tt n 2n a 

T' ^2 - -^ 0: 



3 - 
^3 

0, 

3 - 



27r 
5 ' 
7r 
5' 



0A 



27r 
' 5 ' 



0A = -T 



5' 
' 5 ' 



0,-- 

04 



5' 

— 
5' 
2n 
5 ■ 



Then, M = M4 ■ M3 ■ M2 ■ Mi ■ Mq, and, for the four values of p, 

(1 + cos(p))^ 



Tr(M) 



1 = 1 



2cos(|f), 



so the only possible real value for cos(p) is 

cos(p) = (16 + 16cos(|f ))i/^ - 1 = 2cos2/^(§) - 1. 
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(74) 



(75) 



(76) 



(77) 



(78) 



When q is even, the previous ideas are completely valid. As we saw in the 
previous section, the case p = 1, g = 2 is trivial, so let us consider g = 4. 
Then, we have the following angles: 



(79) 



Ml and M3 are identity matrices, so M = M2 • Mq. Then, for both values 
ofp. 



= 1 : 


n TT n TT 
^0 — -4, ^2 — 4, 


= 3: 


TT n TT 

(70 — 4, (72 — -4- 



Tr(M) = (1 + cos(p))^ - 1 = 1 + 2 cos 



•27^^ 



(80) 



In this case, cos(p) can take two real values, cos(p) = ±(2 + 2 cos(||))"'^/^ — 1, 
but, since p is real, | cos(p)| < 1, so the only valid value for cos(p) is 



cos(p) = (2 + 2cos(|f))i/^ 
If we now choose q = Q, then 



l = 2cos(^)-l. 



p=l: 


n TT n TT 
(71 — -g, ^3 — 2' 


0, 


p = 5 : 


n TT n TT 

01 = g, ^3 = -2' 


0, 



(81) 



(82) 



This angles are exactly the same as in (71), but in a different order. Now 



Mq, M2 and M4 are identity matrices; M = M5 ■ M3 ■ Mi, and the only 



possible real value for cos(p) is again (73). 
Finally, if we choose g = 8, then 



p = l 


^0 = -f , 


O2 = 0, 


^4 = f 


e 


p = 3 


^0 = X' 


O2 = 0, 


^4 = -f 


e 


p = 5 


(70 — X' 


^2=0, 


^'4 = -f 


9 


p = 7 


^0 = -f , 


O2 = 0, 


Q 37r 
(74 — X 


9 



0, 

0, 
0, 
0. 



(83) 



Ml, M3, M5 and M7 are identity matrices, so M = Mg ■ M4 ■ M2 ■ Mq. 
Then, for all the possible values of p we have 



Tr(M) 



;i + cos(p))- 
4 



-1 



l + 2cos(|f). 



and the only real value that satisfies | cos(p)| < 1 is: 

cos(p) = (8 + 8cos(||))V4 - 1 = 2cosV^(f ) - 1. 



(84) 



(85) 
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The previous results strongly suggest that, in general, for any q and for any 
p coprime with it, the only possible real value for cos(p), | cos(p)| < 1, is 

f2cos2/9(^)-l, if g is odd, , , 

cos(p) = { ..Y'[ ' ^ ' (86) 

1 2cos /'^(^) — 1, if g is even; 

remark that if q is odd, the angle p between two adjacent sides is the same 



for q and for 2g. Although giving a universal proof that (86) holds for any 
q goes beyond the scope of this paper, we have checked it for a few more q. 
Moreover, it is absolutely in agreement with our numerical simulations, as 
we will see in Section [5j Therefore, we assume it is valid Vg. 

Once we have found the correct choice of p, we get immediately the value 



for ?/)(0,tpg). Indeed, bearing in mind (59) 



MyTq 



,^ arccos (2 cos^/^l^) - l) , if g is odd. 
^^(0,^= <J A,/| (87) 

-^^ arccos (2 cos '^(-^) — 1) , if g is even. 

It is straightforward to check that '?/'(0, t\^2) = 1- It is also interesting to note 
that 

limV^(0,t,,) = ^(-ln(cos(-)))V^ ^^^^ 

lim ij{0,tpg) = 1. 

Moreover, T, ei and e2 can be completely determined for any p and g, too, 
up to a rigid movement. The same is valid for X, obtained from integrating 
T once. In general, for an arbitrary polygon, it can be pretty complicated 
to specify the correct rigid movement. However, since our initial data is a 
regular planar polygon, the symmetries of the solution are very advantageous. 
The easiest way to understand and to implement numerically the correct 
rotation (although not necessarily the simplest option for symbolic manipu- 
lation) is to work on X, which depends only on T. To that aim, given a time 
tpq, we compute first the associated rotation matrices M^. Then, by means 



of (61), we obtain up to a rotation the piecewise constant vectors T, ei and 



62, which we denote T, ei and 62: 

(89) 





r 


) 


_ ~ (2-Kk-2TT + \T 

~ *^H Mq > 


\ 


e2(i^- 


T 


\ ~ (2TTk-2-K + \T 

\ 621 Mq ) 


) 
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and the case k = can be assigned any value, for instance the identity matrix, 
(T(0~)|ei(0~)|e2(0~))^ = I. Let us underhne that it is not necessary at all 



to consider separately q odd and q even. Indeed, in the case of q even, (63) 



X(^l^) = X(ff ) + i-t(llf-), (90) 



and (65) safely allow us to regard the Mg/2-sided polygon as a degenerated 
Mg-sided polygon. This is extremely useful for implementation purposes. 

Once T has been obtained, X, which is X up to a rotation, is computed 
recursively: 

^^- , -.^, 

Mq ' ^Mq' ' Mq ^ Mq 

where, again, X(0) can be given any value, for instance X(0) = (0, 0, 0)"^. 
To obtain the correct rotation for X, we use the symmetries of the polygon, 
as explained in Section |3j In particular, we use the facts that, for any time 
t, X.[27rk/M), k = 0, . . . , M — 1, have to be coplanar and lay on a plane 
orthogonal to the 2;-axis; and that X(27r/M) — X(0) is a positive multiple of 
(1, 0, 0)^. This can be done efficiently as follows: 

1 Tomniite v+ = X(2VAf)-X(0) ^ X(0)-X(-2^/M) 

^ ||X(27r/M)-X(0)|| ' ||X(0)-X(-27r/M)ir 

2. Compute w = v^ A v+. 

3. Compute the scalar product between w and (0, 0, 1)"^, w ■ (0, 0, 1)-^ = 

4. If W3 = 0, Ri is the identity matrix. If not, Ri is the rotation matrix 
that induces a rotation of aiccoslw^) degrees around the axis given by 



the vector 



wA(0,0,l)-' 



||wA(0,0,l)^||- 

5. Compute v+g^ = Ri ■ v+. 

6. Compute the rotation matrix R2 that induces a rotation of arccos(v+g^- 
(1,0,0)^) degrees around the z-axis. 

7. Compute the sought rotation matrix, R = R2 ■ Ri. 

8. Update T = R ■ t and X = R • X. 

By means of the previous algorithm, we have computed the correctly rotated 
values of X and T. Additionally, for X, we still have to specify a translation. 
Again, bearing in mind the symmetries of the problem, we translate X in such 
a way that its mass center, given by the mean of X(27rfc/M), k = 0, . . . , M—1, 
is in the z-axis. In this way, we have completely specified the position of X, 
up to a vertical translation. 
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4. Numerical method 

In this section, we simulate numerically (tsl) and Q taking respectively 



(28) and (29) as initial data, i.e., X(s,0) is a planar regular polygon of M 
sides and length 2ti. 

There have been a couple of papers devoted to reproducing numeri- 
cally the self-similar solutions of the Schrodinger map (ll]). In [8], a Crank- 
Nicholson scheme was consider, together with a finite-difference discretization 
of Tss- Later on, in [U], both finite-difference discretizations and pseudo- 
spectral discretizations were considered to simulate ^ and its equivalent on 
the hyperbolic plane. In particular, given the asymptotical structure of the 
self-similar solutions of (|4]), a truncated domain s G [—L^L]^ L 3> 1, with a 
grid based on the Chebyshev nodes was found to be very convenient to rep- 
resent those solutions. Nevertheless, due to the clustering of the Chebyshev 
nodes, an explicit scheme to advance in time implied the undesirable restric- 
tion |At| = 0{1/N^), where A^ is the number of nodes. In order to solve that, 
it was chosen to work with the stereographic projection of T = (Ti,T2,T^) 
over C, 

/ , f Ti T2 \ 

z = X + iy = (x,y) = — -, — - , (91) 

which transforms (111) into a nonlinear Schrodinger equation: 

zt = izss - I .^ zl (92) 

"*" I I 

The advantage of this equation, as opposed to (|4]), is that the higher-order 
term Zgs can be treated implicitly, therefore eliminating or at least reducing 
significantly the restrictions on At. 



Working with (92) was a very adequate choice for the purposes of [9]. 
There, a second-order semi-implicit backward differentiation formula was 
chosen. This scheme was very stable because it imposed a very strong decay 
in the high frequency modes; moreover, its low order was easily and effec- 
tively compensated by the use of an adaptive method, both in space and in 
time. 



However, in this paper, working with (92) does not seem to be such a 
good option. Indeed, unlike in j9j, where an extremely high accuracy was 
required only for < t <^ 1, we are now interested in the behaviour of T 
at all times, for which a second-order scheme seems a very poor option. 
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After many unsuccessful attempts, which have greatly delayed the redaction 
of this paper, we admit that we have been unable to find a good higher-order 
scheme for (92). Furthermore, we are also interested in the evolution of X, 
for which working with (92) is of limited help. Bearing in mind the previous 
arguments, we have opted to work directly with ^ or, more precisely, with 
a combination of ^ and (|4]): 

-^^^- (93) 

= TAT„. ^ ^ 



In order to simulate numerically ( 93 ) , we have chosen a fourth-order Runge- 
Kutta scheme in time: 

Ax = T(") a T("), At = T(") A T^^), T^^^ = T^") + f At, 
Bx = T(^) a T(^), Bt = T(^) a Tif, T(^) = T(") + f Bt, 
Cx = T(^) a Tp), Ct = T(^) a T^f ), T(^) = T(") + AtCr, 
D;, = T(^)ATf), DT = T(^)ATi^), 
X{"+i) = x(") + f (Ax + 2Bx + 2Cx + Dx), 

X = t(") + ^(At + 2Bt + 2Ct + Dt), t("+^) = ^. (94) 

The last line guarantees that T*^"^ G S^, Vn. Additionally, we can also project 
T^i^)^ T^^\ T^'"' onto the unit sphere, but we have noticed no significant 
improvement in our numerical results. 

We have combined the previous Runge-Kutta scheme with a pseudo- 
spectral discretization directly in space. More precisely, since we deal with 
periodic solutions in s G [0,27r), we have simulated the evolution of X = 
(Xi,X2,X3) and T = (Ti,T2,T3) at N equally spaced nodes Sj = 2nj/N, 
j = 0, . . . ,N — 1. In order to compute T^, T^^, etc., we remember that, 
given a periodical function /(s) evaluated at Sj, its derivatives at Sj can be 
spectrally approximated as 

N/2-1 N/2-1 

fs{s,)= Yl ^fc/>)e'"^'/'^, fss{s,) = - J2 ^V>)e''^'^''/^, (95) 

k=-N/2 k=-N/2 

where 

7V-1 
j=0 
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Both (95) and (96), wich are respectively inverse and direct discrete Fourier 



transforms (DFT), can be computed efficiently by means of the Fast Fourier 
Transform (FFT) algorithm [2l]. Moreover, since X and T are invariant by 



rotations of 27r/M degrees around the z-axis, as shown in (32), we can reduce 



the computational cost of (95) and (96) to a DFT of N/M elements, instead 



of a DFT of A^ elements. In what follows, we explain how to apply this idea 



to T. This is valid, with no change, for any vector satisfying (32), i.e., X, 
T„T 



ss, etc. First, let us consider the ffist two components of T. Denoting 



Ti{sj,t) + iT2{sj, t), we have 

N-l 

Z{k) = Y, ^(s,-)e-'"^''^/^ 

j=0 

A/-1 N/M-1 

1=0 j=0 
'M-l \ N/M-1 



=0 / i=0 

N/M-1 

M J2 Z{sj)e-'^^'^^l^ , if fc = 1 mod M, 

j=o 
0, if A; ^ 1 mod M. 



(97) 



Therefore, the only non-zero Zilz) are those with k = 1 mod M. Moreover, 
they can be obtained by a DFT of N/M elements instead of a DFT of A^ 
elements: 

N/M-1 

Z{Mk + 1) = M ^ Z(s,-)e-2"*^'(^^*^+')/^ 
i=o 

N/M-1 

= M J2 [e-2"*^'/^Z(s,)] e-2-^'*=/W^). (98) 

i=o 
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Identically for the third component of T, 

Af-l 



7—0 


^-2nijk/N 






M-l N/M-1 
/=0 i=0 


Tz{Sj+i{^]^/M))e 


2m{j+l(N/M))k/N 




/M-l \ N/M-1 
\ 1=0 J j=0 


^,.y-2.^jk/N 




N/M- 


L 


if A; = mo 


dM, 


3-0 

0, 




if fc ^ mo 


d M. 



(99) 



Therefore, the only non-zero T-i{k) are those with /c = mod M and, again, 
they can be obtained by a DFT of N/M elements instead of a DFT of N 
elements: 

N/M-1 

%{Mk) = M Y^ T3(sj)e-2"'^'^/W*^). (100) 

3=0 

Furthermore, since T^{s) is real, this last expression can be further simplified 
to a DFT of N/{2M) elements. 



Combining (94) with a pseudo-spectral discretization in space, we obtain 



a scheme with a time-step restriction on At that appears to be of the form 
At < C/N^, with C ~ 11.3. However, for a fixed M, we only want to 
simulate X and T for t G [0,2n/M^]. Denoting Atmax the biggest At that 
makes ( [94| ) stable, the smallest number of time-steps A^j needed to reach 
t = 27r/M2 is of at least 



N, > 



"27r/M2" 

^l^max 


= 


"27r/M2" 

c/m 



N\- 

m) 



27r' 
~C 



(101) 



From this equation, together with (98) and (100), it follows that, for a given 



r G N, if we take A^ = 2^M, i.e., the computational cost for simulating (94) 



is exactly the same for any number of initial sides M, no matter how big it 
is. This important fact enables us to make consistent comparisons between 
different M, because we are considering an equivalent spatial resolution, i.e. 
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we are simulating the same number of points N/M. We take N/M a power of 
two to take the greatest possible advantage of the FFT. Obviously, for N/M 
constant, At decreases as (9(1/M^); therefore, the most accurate results will 
be expected for the largest M , as we will see in our numerical experiments. 
In the following section, we simulate X and T for initial data with dif- 



ferent M. As we will see, despite its simplicity, (94) gives surprisingly good 
results. 

5. Numerical experiments 

In Section |3} in order to construct algebraically the evolution of a regular 
polygon at rational times, we have done some very strong assumptions, the 
most important one being uniqueness. However, as we will see in this section, 
the numerical experiments are in complete agreement with their theoretical 
predictions. To illustrate this, we have simulated the evolution of X and T 



by means of (94), taking regular polygons with different numbers of sides M 



as initial data, and making X and T evolve until t = 2tt/M'^. The initial 



data X(s, 0) and T(s,0) are given respectively by (28) and (29); obviously, 
in the case of X, the non- vertex points are immediately computed by linear 
interpolation. We have taken N equally-spaced nodes Sj = 2tt j/N G [0,27r), 
0, . . . , A^ — 1; nevertheless, bearing in mind the symmetries of X and T as 
described in the previous section, we only have to describe the evolution of 
X and T at the first N/M nodes Sj G [0, 2tt/M), j = 0, . . . , N/M - 1. We 
have divided the time-interval [0, 27r/M^] in Nt equally-spaced time steps of 
length At = -j^r^q-- Therefore, during the simulation of (94), we obtain X^"') 



and T^'^^ at t^"-* = riAt, ra = 1, . . . , Nt. In our experiments, bearing in mind 



( [101| ), we have chosen Nt = 151200 ■ 4^ for N/M = 512 • 2\ 

Remember that, in Section |3} we had constructed algebraically X up to 
a vertical movement. Therefore, in order to completely specify X at a given 
time, we would need to give the height of one point or, more conveniently, 
the height h{t) of the mass center, which is precisely the mean of all the 
values X2,{sj,t): 

1 ""-' m""^'^ 

M^) = ]^E^3(^„t) = - Y. Ms„t). (102) 

i=o i=o 
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Figure 2: Evolution of h'{t) = niean(X3^f) against t, for an initial triangle, M = 3, 
iV = 3 X 4096 nodes. 



On the one hand, h'{t), computed numerically by 

N-l 
j=0 

- 5^(Xi,,(s„t)X2,,,(s,-,t) -X2,,(s,-,t)Xi,,,(s,-,t)), (103) 



N 

j=0 

is far from constant and has a very singular shape, as shown in Figure |2} 
where we have plotted the numerical evolution of h'{t), at all t'-"^ = nAt, for 
an initial triangle, M = 3, A^ = 3 x 4096 nodes. 

On the other hand, when integrating h'{t), the oscillations cancel and 
completely disappear. Furthermore, h(t) can be approximated with great 
accuracy by means of a constant multiplied by t; more precisely. 



hit) . ^^t = CMt, 



(104) 



where cm = h{^)/ ^ is the mean speed. In Table ^ we give the maximum 
discrepancy between h{t) and its linear approximation cm^, i-e., max„ \h{t^'^^) — 
Chit I • We have considered different M and different numbers of nodes N] 
the errors are very small and they seem to decrease as 0{1/N). This gives 
very strong evidence that h{t) is linear or, at least, quasi-linear. Therefore, 
in practice, we can safely assume that h{t) = cut- Observe that cm, which 
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M 


N/M = 512 


N/M = 1024 


N/M = 2048 


N/M = 4096 


Cm 


3 


4.3096 


10-5 


2.1206 


10-5 


1.0886 


10-5 


5.7953 


10-6 


0.7645 


4 


1.2398 


10-6 


6.1344 


10-6 


3.2140 


10-6 


1.7316 


10-6 


0.8827 


5 


4.8504 


10-6 


2.4191 


10-6 


1.2807 


10-6 


6.9338 


10-^ 


0.9286 


6 


2.2848 


10-6 


1.1441 


10-6 


6.0905 


10-^ 


3.3044 


10-^ 


0.9517 


7 


1.2167 


10-6 


6.1060 


10-^ 


3.2607 


10-^ 


1.7710 


10-^ 


0.9651 


8 


7.0721 


10-^ 


3.5594 


10-'^ 


1.9014 


10-^ 


1.0333 


10-'^ 


0.9735 


9 


4.3905 


10-^ 


2.2140 


10-'^ 


1.1828 


10-^ 


6.4303 


io-« 


0.9792 


10 


2.8697 


10-'^ 


1.4489 


10-^ 


7.7407 


io-« 


4.2093 


io-« 


0.9832 



Table 1: max„ |ft,(i'"'') — cm^'-"''!, for different numbers of initial sides M, and of nodes N. 
^(n) = nAt; Cm = ^{'m^)/ 'm^- The values corresponding to cm have been calculated with 
N/M = 4096. 



is also offered in Table [I] for N/M = 4096, grows with M, tending to 1, as 
M tends to infinity, i.e. as X(s,0) tends to a circle. 

In order to compare the values of X obtained numerically, which we la- 
bel X„um, with those obtained algebraically, which we label ^aig, we have 
subtracted the vertical position of the center of mass from X„„m, i.e., we 
have analyzed the agreement between X^^^ — cm^(0, 0, 1)^ and X^/g. More 



precisely, we have computed maxm(maxj ||X 
Xazg(s„t(™))||), where s, = 2ttj/N, 
Euclidean distance; and t^™'^ 







t^™)) -CAft(") (0,0,1)^ - 
■ II denotes the 
Notice that, in 



,N/M-1; 
1260. 



, J 

2tt m 

., ^^^^ V — j^p ■ 1260' 

Sectional we have constructed Xa/g(s,t) only at the vertices, so the non- 
vertex points ^aig{sj,t) are computed again by linear interpolation. Observe 
also that comparing ^num and X^jg at A^^ + 1 time-instants would have been 
computationally unrealistic. Instead, we have chosen 1260 + 1 time-instants, 
because 1260 is a number large enough for our purposes, and with a conve- 
nient factorization, 1260 = 2^ ■ 3^ ■ 5 ■ 7. 

The maximum value of ||X„„„(sj, t^™)) -CAft("')(0, 0, 1)^-X„ig(sj, t(™))|| 
is given in Table [21 for different N/M and M. Bearing in mind that we are 
comparing X„„m and '^aig globally for a large number of nodes and time- 
instants, and that max ||Xa;g|| > 1, VM, the results are, in our opinion, very 
remarkable, and strongly suggest that there is convergence, or at least an 
extremely good agreement, between both approaches, as A^ — )■ oo; moreover, 
since X^^^ is periodical in time, with time-period 27r/M^, they also suggest 
that X is periodical, or at least quasi-periodical, in time with that period, 
up to a vertical movement. All this is also supported by the fact that the 
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M 


NIM = 512 


NIM = 1024 


N/M = 2048 


N/M = 4096 


3 


2.4847 


10-3 


1.3841 


10-3 


8.1211 


10-4 


4.9718 


10-4 


4 


1.1221 


10-3 


6.9665 


10-4 


4.2717 


10-4 


2.5917 


10-4 


5 


6.8414 


10-^ 


4.2545 


10-4 


2.6125 


10-4 


1.5874 


10-4 


6 


4.6057 


10-4 


2.8717 


10-4 


1.7670 


10-4 


1.0754 


10-4 


7 


3.3170 


10-4 


2.0724 


10-4 


1.2772 


10-4 


7.7832 


10-5 


8 


2.5059 


10-4 


1.5680 


10-4 


9.6744 


10-5 


5.9010 


10-5 


9 


1.9616 


10-4 


1.2288 


10-4 


7.5878 


10-5 


4.6313 


10-5 


10 


1.5782 


10-4 


9.8943 


10-5 


6.1137 


10-5 


3.7334 


10-5 



Table 2: max™(niax, ||X„„^(s,, i^"') - CMt("')(0,0, 1)^ - X,,g(s,,i("))|i), where s^ 



^, J = 0,---,7V/M 
m = 0,...,1260. 



1 



denotes the Euchdean distance; and t*^™) 



2tt 



1260' 



errors decreases as M increases. Indeed, from (101), we are taking the same 
number of time-steps Nf for each M, but At = jfljj-, i.e.. At = C(l/M^), 
which explains those smaller values for the last rows of Table [2] and also of 
Table [H 

On the other hand, comparing the algebraically constructed Taig with 
the numerically obtained Tnum is more problematic. For example, in Figure 
|3| we have compared Tnum (left) with Taig (right), for an initial triangle, 
M = 3, at ti,3 = If, N/M = 4096. The exact value of T at that time is 
given by Taig, which, by construction, is a piecewise constant function with 



exactly qM = 9 pieces, that we denote Tj, i = 1, 
values, in this case, have an easy algebraic expression: 



, 9, and whose explicit 



^-1 



^-1 



Ti 



/4- 1 
1-^ 




(105) 



T4, T5 and Tg are respectively Ti, T2 and T3 rotated ^ degrees around 
the 2;-axis; T7, Tg and Tg are respectively Ti, T2 and T3 rotated ^ degrees 



■7, J-8 

around the 2;-axis. 



Observe that T^ 



exhibits a Gibbs-type phenomenon. However, in Tnum, if we take the central 



is almost identical to Taig, but clearly 
)wever, in Tnum, if we take the central 

values of each interval (altogether 27 circles, indicated in Figure |2] with a 

black circle) 

maximum error equal to 6.3869 



and compare them componentwise with (105) 

9 



we obtain a 
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Figure 3: Tnum versus TaZg, for M = 3, at ii_3 = |f • T'l appears in blue, T2 in green, T3 
in red. In Tnmn, the Gibbs phenomenon is clearly visible. The black circles denote the 
points chosen for the comparisons. 

5.1. X(0,t) and Riemann's non-differentiahle function 

Bearing in mind the mirror symmetries of the initial data, as illustrated 
for the pentagon in FigureOl we conclude that X(7rfc/M, t), fc = 0, . . . , 2M— 1, 
lives in a plane that contains the z-axis. Indeed, s = nk/M are the only 
points such that X(s, t) describe a planar curve. In what follows, we will 
describe X(0,t), although all said here is immediately applicable to any s = 
nk/M. Since X(0,t) = (Xi(0, t),X2(0, t),X3(0,t)) is planar, bearing in mind 
that Xi(0, t) < 0, X2(0, t) < 0, we rotate X(0, t) clockwise it/2-7i/M degrees 
around the z-axis, until it lays on the plane OYZ, which we identify with C. 
Then X(0, t) becomes 



zit) = -||(Xi(0,t),X2(0,t))|| +zX3(0,t). 



(106) 



In Figure El we have plotted the numerically obtained z{t) for different M. 
Besides the conspicuous fractal character of the curves, which immediately 
reminds us of the works pil| [TT] , their most striking feature (see Figure [s]) 
is how much z{t) — icMt, i.e., z(t) without the vertical movement, resembles 
the graph of 



m = E 



°0 gTTifc^i 



fc=l 



ink"^ 



te[0,2]. 



(107) 



(f)(t) was conveniently used in ^/Tj, in order to study Riemann's non-differen- 
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Figure 4: z{t), for M = 3, 4, 5, 6, 7, 8. N/M = 4096. We have plotted the z{t) correspond- 
ing to aU the Nt + 1 points X*^"' (0) obtained during the execution of ( 94 ) . 



tiable function, which is precisely the real part of 0(t), i.e., 

sin(7r/c^t) 



m = E 



k=l 



irk? 



(108) 



f{t) is non-differentiable everywhere, except at those rational points t = p/q, 
with p and q both odd. On the other hand, when constructing algebraically 
X(0,t) for a given t, we have observed that the only times at which X(0,t) 
has no corner are of the form tpq = jjf-, with p odd and q = 2 mod 4. 
Therefore, in order to compare (j){t) and X(0,t), we have to redefine slightly 



0(t) in (107): 



0(t) 



TT 

6 



-E 



°° ^2nik^t 



k=l 



TikP- 



tG[0,l], 



(109) 



i.e., we multiply (p in (107) by — z, in order to orientate it correctly; then 



we translate it, so 0(0) = 0(1) = 0, and, finally, we change its period from 

t e [0,2] tot G [0,1]. 



On the other hand, given M, we compare in (109) with 



ZM{t) = z{^)-z{Q)-lCM^2, te[0,l]. 



;iio) 
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Figure 5: (/)(i) (red), as defined in (1091, against ZM{t) (blue), for M ~ 3,4,5,6,7,8. To 
compute (j){t), we took J2k=o^ ■ "^^ compute, ZM{t), we took N/AI = 4096. We have 
plotted the ZM{t) corresponding to all the Nt + 1 points X(")(0) obtained during the 
execution of (94). 



In Figure |5| we have plotted 4'{t) (red), as defined in (109), versus zm(^), 
for different M. Althougli tfie figures are not identical, they are extremely 
similar. Moreover, if we scale them accordingly, it seems that the bigger M 
grows, the more similar to 0(t) is ZMit). In what follows, we will give some 
evidence that this indeed happens. 

Since we have obtained numerically -2a/ (t) at A^i + 1 points t G [0,1], 
we have to evaluate 0(t) precisely at those points. Then, we scale ZM{t), 
approximating the scale factor in this way: 

m.axRe((h(t)) 

, te[o,i] 



\M 



maxReizMit)) 

te[o,i] 
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M 


4>{t) - XMZM{t)\ L°° 


3 


9.2128 


10-^ 


4 


4.7210 


10-^ 


5 


2.9124 


10-^ 


6 


1.9899 


10-^ 


7 


1.4517 


10-^ 


8 


1.1091 


10-^ 


9 


8.7713 


10"=^ 


10 


7.1269 


10-=^ 



Table 3: Maximum error between 0(i) and ZMit) (scaled). To compute ()){t), we took 
J2k=o~ ■ "^^ compute ZM{t), we took N/M — 4096. We have considered the ZM{t) 
corresponding to all the Nt + I points X("^(0) obtained during the execution of (94). We 
have computed algebraically 4>{t) at those times. 



Finally, we compute ||0(t) — Am^m(^)||l°°- Table |3] shows the errors obtained 
in this way for several M. Bearing in mind that maxtg[o,i] \4'{t)\ ~ 0.9187, 
that the scahng Xm was not intended to be optimal, and that we are compar- 
ing 9676800 + 1 different points one by one, the results give a pretty strong 
evidence of that convergence. This gives strong support to the idea that zm 
is also a multifractal, VM > 3. 

As we mentioned in the introduction, the proof that f{t) is a multifractal 
was done in [T8j. In this respect, it is fundamental the identity given for 0, 
as defined in (|107), which states that 



0(t) = 0(t 



pg. 



+ e"^"/^g-i/2(^ 



tpqY'^ + lower-order terms. 



;iii) 



with p, g G Z, g > 0, gcd(p, g) = 1, 



771 



7n{tpq) G Z/8Z. This identity is 



proved in Theorem 4.2 in [17j. Of particular relevance is the fact that the 
Holder exponent is 1/2, i.e., that |0(t)— 0(tpg)| = q''^^'^\t—tpq\^/'^+ lower-order 
terms. 

In order to prove analytically that ZM{t), is a multifractal, we would need 
We have made some numerical experiments (see for 
3, p = 1, g = 5), and all of them give strong 



an equivalent of (111) 

example Figure l6| for M 

evidence that the Holder exponent of z{t) is 1/2 for rational times, i.e., that 

\ZMit) - ZMitpq)\ =Oi\t-tpq\^/^). 

Nevertheless, how the constants depend of the denominator g, something 
which is a fundamental ingredient in the arguments in [TTj and in 
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Figure 6: |z(t) — z(ti.5)| versus (i — ii.5)^/^, for M = 3. The asymptotically linear relation 
between both quantities as i — >■ ii^s is evident. 

unclear. This question deserves a much more detailed analysis that we plan 
to make in a forthcoming paper. 



5.2. T{s,tpg),forq:^l 

In the previous subsection, we have given evidence of the multifractal 
character of X(0,t). However, in [9], fractal-like phenomena were observed 
also on T, when imposing periodic boundary conditions. In the following 
lines, we will show that similar observations are valid for periodic boundary 
conditions, too. 

As we have seen in this paper, X can be recovered algebraically up to 
a vertical movement, while T can be completely recovered algebraically at 
rational times. Furthermore, Taig is really a piecewise constant function, with 
no noise associated to the Gibbs phenomenon. This gives us a very powerful 
tool to progress in the understanding of T (and also of X), avoiding numerical 
simulations at all. 

A very interesting question is, given a rational time tpg = j^, with q 
small, corresponding in X to a polygon of Mq or Mq/2 sides, and in T to a 
piecewise continuous function with Mq or Mq/2 jumps, what happens at a 
time t + e, \e\ <^ 1? 

Let us take e = -^^^ 

gcd(g, q') = 1, then, 2 + ^ 
of qq' or qq' /2 sides 



with g' ^ 1, in order that \e\ ^ 1. Assuming that 
1 _ 22_±i gQ ^ _|_ ^ would correspond to a polygon 



qq' 



In Figure we have have plotted ^aig and Taig, for M 



3, at t — -g-(4 + 
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Figure 7: X^,,, and T^/g, at i = ^(i + j^) 



1 N _ 27r 50003 



While visually there is no difference whatsoever between 



Xa/„ at ti 4 = if ■ i and X^,, at t = ^f ■ W^, we do not have a skew polygon 



49999 ^ 9 199996 ' 

•-alg at ti4^ ^ "g" ■ 4 anu -^^alg ai fc ^ -g- ■ Jggggg; 

with 6 sides, but a skew polygon with Mq/2 = 299994 sides that closely 
resembles a polygon with 6 sides. On the other hand, the corresponding Taig 
are very different. Indeed, in the plot of Taig, we can clearly appreciate six 
spiral-like structures whose centers are precisely the six constant values of T 
at ti_4. Nevertheless, it is important to underline that these structures are 
not really spirals, but the plot of the 299994 different values taking by T, that 
closely resembles a curve with six spirals. Furthermore, these spirals remind 
us of the Cornu spirals that appeared in [6]. An open question that arises 
naturally is up to what extent the multiple-corner problem can be explained 
as a sum of several one-corner problems. 

Another interesting question is what happens if we take a time tpg with 
large q, such that there is no pair {p',q'), with both q' and |- — ^| small. In 
this case, the situation is very different. In Figure |8l we have plotted ^aig 
and Taig for M = 3, at t = f (i + ^ + ^) = f ■ |||. While the left-hand 
side is not so different from the left-hand side of Figure [7| the right-hand side 
is a set of Mq/2 = 98646 points in §^ that creates a spectacular fractality 
sensation with spiral-like structures at three or four different scales. This can 
be better appreciated in Figure |9| where we have plotted the stereographic 
projection of T onto C 

In general, a time tpg, with g ^ 1, can be regarded as an approximation 
of an irrational time. Trying to fully understand the phenomena exhibited 
by T(s,t) at those times is a challenging question that we postpone for the 
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Figure 8: Xaig and Taig, at t 
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Figure 9: Stereographic projection of the right-hand side of Figure [8) 



future. 
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6. Conclusions 

In this paper, we liave studied the evolution of ^ and Q, for a regular 
planar polygonal of M sides as initial datum. The algebraic calculations, 
backed by complete numerical simulations, suggest very strongly that X(s, t) 
is a polygon at times which are rational multiples of 27r/M^, i.e., tpq = 
(27r/M^)(p/g), with the number of sides depending on q, while T{s,t) is 
piecewise constant at those times. 

We have obtained a striking connection between X(0, t) and the so-called 
Riemann's non-differentiable function. In [TS], S. Jaffard proved that this 
function is an example of multifractal whose spectrum of singularities satisfies 
the Frisch-Parisi conjecture. Although there is strong numerical evidence 
that X(0, t) is also a multifractal, an analytical proof seems to be challenging. 
In fact, a first step in this direction is giving sense to our solutions from an 
analytical point of view. Moreover, while our X(s,t), obtained algebraically 
for rational times, can be extended by continuity to all t G M, it is not clear 
how to give sense to T(s,t) at irrational times. 

In the future, we also plan to extend these ideas to arbitrary polygons, 
and to do a more detailed study on fractality. For such purposes, a more 
complete algebraically description of X and T is no doubt required. More 



precisely, we would aim at a level of detail similar to that in (86). 
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Appendix A. Generalized Quadratic Gaufi Sums 

The generalized quadratic GauB sums are defined by 

|c|-l 

G(a, 6, c) = ^ e2"'('^''+^')/^ (A.l) 

1=0 

for given integers a, b, c, with c 7^ 0. From now on, we assume c > 0, and 
gcd(a, c) = 1, which are the cases dealt with in this paper. 
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The value of these sums was calculated for the first time by GauB for the 
case 6 = [25]. More precisely, given two integers a, c such that gcd(a, c) = 1, 
then 



c-l 



G(a,0,c) = ^, 



^2mafi /c 



'(f)(l + i'^)yc, ifc = mod 4, 

(f)v/^, ifc^l mod 4, 

0, if c = 2 mod 4, 

X-)iy/c, ifc = 3 mod 4, 



(A.2) 



where (-) denotes the Jacobi symbol. 

As we will see in the following lines, the generalized quadratic GauB sums 
can be reduced to normal quadratic GauB sums by completing the square. 
Since they are multiplicative, i.e.. 



G(a, 6, cd) = G{ac, b, d)G{ad, b, c), with gcd(c, d) 



(A.3) 



we can assume without loss of generality that c is either odd or a power of 
two. If c is odd, let us find a certain ip{a) such that 4aip{a) = 1 mod c. 
This '?/'(a) is unique modulo c and is precisely the inverse of 4a in Z/cZ; since 
gcd(4a,c) = 1, its existence is guaranteed by Bezout's lemma, and it can 
be efficiently computed, for instance, by the extended Euclidean algorithm. 
Then, we have 



c-l 



G{a, b, c) = y^ e 



2TTi{al^+4:a^p(a)bl+i>{a)b^-tp{a)b^)/c 



1=0 



-2mip{a)b'^ / c 



c-l 

E 

c-l 



^27rJ(ai2+4a?/)(a)W+4ai/i(a)V'(a)fe^)/c 



-2-nitp(a)b'^ / c ST^ 2nia{l+2i}{a)bf/c 

1=0 
c-l 
-2niil>{a)b'^/c Sr^ 2mal^/c 



■ ST^ 2-Kia 



1=0 

-2-Kiil){a)b'^ / c/ a 

\ c 

-2-Ki'il){a)bi^ /cf a 



)a/c, if c = 1 mod 4, 
)iy/c, if c = 3 mod 4. 



(A.4) 
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If c is a power of two, we first observe that 

c/2-l c/2-1 

G{a, b,c) = y^ f^2niial^+bl)/c _^ V^ g27ri(a(i+c/2)2+6(i+c/2))/c 

1=0 1=0 

c/2-1 
= \l + (_l)«(c/2)+61 y^ g27ri(ai2+bO/c 

1=0 

'O, ifa(c/2)^6 mod 2, 

2 E^fo"^ e^^i(ai^+bi)/c^ if ^(^/2) = 6 mod 2; 

therefore, we have two cases: c = 2 and c > 2. If c = 2, then a is odd, so 

^/ r ON 1 0' if ^ is even, 

G(a,6,2) = <^ ' (A.5) 

12, it IS odd. 

On the other hand, if c > 2, G{a, b, c) = 0, if 6 is odd. Let us suppose b is 
even. Then, since a is odd, we follow the previous reasoning to find a certain 
ip{a) such that aip{a) = 1 mod c, so 

c-l 

2Tri{al'^+atl>{a)2{b/2)l+i;{a){b/2f-ip{a){b/2)'^)/c 



G{a, b, c) = y ^e 

E 

c-l 



z=o 

c-l 

TTJ^Wb^/CSc) V^ 2m(al^+2aij(a)(b/2)l+a'tl){a)i}{a){b/2)'^)/c 

1=0 
c-l 
g 7r#(a)b2/(2c) Y^ g27ria(«+V(a)(fe/2))Vc 

/=0 
c-l 



-TTJVWb^/CSc) V^ 2TTiafi/c 

1=0 

-nii,{a)by{2c)^cj^^^-a^^_ (A.6) 



In general, given an arbitrary c G N, we will factorize it as c = 2'^c' and use 
the multiplicative character of the generalized quadratic Gaussian sums to 
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calculate G{a, b, c). For instance, if c is even, but c/2 is odd, then 

G{a,b,c) = G{a,b,2{c/2)) 

= G{a{c/2),b,2)G{2a,b,c/2) 

0, if a(c/2) ^b mod 2, , , „, 

V / 7 7- (A.7) 

2^(20,6,0/2), ifa(c/2) = 6 mod 2, 

i.e, G{a, b, c) = 0, if 6 is even, etc. 

To conclude this appendix, we will mention the explicit value of \G{a, b,c)\, 
deduced from the previous calculations, which is of especial relevance in this 
paper. If gcd(a, c) = 1, then 

{^/c, if c is odd, 
y/2c, if c is even and c/2 = b mod 2, (A. 8) 

0, if c is even and c/2 ^ 6 mod 2. 
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